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We have derived and implemented a stress tensor formulation for the van der Waals density 
functional (vdW-DF) with spin-polarization-dependent gradient correction (GC) recently 
proposed by the authors [J. Phys. Soc. Jpn. 82, 093701 (2013)] and applied it to non¬ 
magnetic and magnetic molecular crystals under ambient condition. We found that the cell 
parameters of the molecular crystals obtained with vdW-DF show an overall improvement 
compared with those obtained using local density and generalized gradient approximations. 

In particular, the original vdW-DF with GC gives the equilibrium structural parameters of 
solid oxygen in the a-phase, which are in good agreement with the experiment. 
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I. INTRODUCTION 

The oxygen molecule has a spin-triplet ground state, and the competing magnetic and molecular 
interactions result in a rich variety of structural, electronic, and magnetic phases under wide 
temperature and pressure ranges.— Much effort has been devoted to clarifying the magnetism and 
atomic structure of solid oxygen.— The magnetism of solid oxygen is particularly interesting, as 
the magnetic moment is sensitive to the relative distance between molecules owing the short-range 
feature of the direct magnetic interaction.— Thus, accurate prediction of the crystal structure is 
one of the key issues for understanding the magnetism in solid oxygen. Theoretical works on solid 
oxygen in the literature are mostly for those under a high pressure, presumably because under 
a relatively low pressure, the vdW interaction plays an important role, where the conventional 
semilocal density functional theory (DFT) fails to describe it properly. 

An accurate description of the vdW interaction is, however, challenging when using the elec¬ 
tronic structure method, particularly, DFT within the semilocal approximation, which is known 
to describe the dispersion interaction poorly. In quantum chemistry, there are well-established 
methods, such as the second-order Mpller-Plesset perturbation theory (MP2), coupled-cluster cal¬ 
culations with single, double, and perturbative triple excitations [CCSD(T)], and the symmetry- 
adopted perturbation theory (SAPT), which can describe the vdW interaction. The quantum 
Monte-Carlo (QMC) method and the adiabatic connection fluctuation dissipation theorem within 
the random phase approximation (ACFDT-RPA) are also known to describe the vdW interaction 
accurately. However, because these highly accurate methods are computationally very demanding 
and applicable only to small systems, practical methods based on DFT are desirable for the simula¬ 
tions of complex systems. In this sense, semi-empirical dispersion-corrected DFT— is widely used, 
and recent developments^^— allow more accurate calculations of the vdW interaction with an 
almost negligible overhead. A more elaborate approach based on the maximally localized Wannier 
function was also proposed to correct for the missing dispersion interaction in DFT.— 

Among other approaches, the van der Waals density functional (vdW-DF), which was developed 
by Rydberg et al— and Dion a/.,— has attracted much attention, because it does not require 
empirical parameters and depends only on charge density and its gradient, thereby allowing one to 
perform large-scale calculation at a modest computational cost. There have been applications of 
vdW-DF to several systems, and to extend the applicability of the methods, there have been several 
proposals for more accurate vdW-DFs.— However, because the vdW interaction is not directly 
related to spin polarization, a truly spin-polarized version of vdW-DF has not yet been developed. 
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except for the one proposed by Vydrov and Van Voorhis.— Thus, the application of vdW-DFs to 
magnetic systems is quite limited. In our previous work,— we proposed a practical approach to a 
magnetic system within the framework of vdW-DF, in which spin-polarization-dependent gradient 
correction is added to the local correlation energy and potential, instead of developing the spin- 
polarized version of the nonlocal correlation. The method has been applied to an oxygen dimer, 
and shown to be applicable to magnetic systems. 

In this work, we have studied solid oxygen at ambient pressure in the a phase (a-02) using 
our vdW-DF for magnetic materials. We implemented the stress tensor arising from the nonlo¬ 
cal correlation to optimize the cell parameters of 0-02 and used several vdW-DFs with different 
exchange and nonlocal correlation. We found that one of the vdW-DFs shows a significant im¬ 
provement on the structural properties over the conventional local density approximation (LDA) 
or generalized gradient approximation (GGA). The remaining part of this paper is organized as 
follows. In Sect. El we describe our implementation of vdW-DF and the computational detail. To 
test our implementation of vdW-DF, we first applied our code to solid argon, graphite, selenium, 
and dry ice (solid carbon dioxide), in which vdW forces play an important role in the binding, and 
the results are given in Sects. IIII AlUTlDI The results for solid oxygen are discussed in detail in 
Sect. IIII El The summary for this paper is given in Sect. IIVI 

II. METHOD 

A. van der Waals density functional 

The exchange-correlation energy in vdW-DF is giveir^ by 

= E, + + Ef, (1) 

where Ex, E'°'^, and E'f^ are the exchange energy, short-range local correlation energy, and nonlocal 
correlation energy, respectively. E°* is expressed as a double integral of the nonlocal interaction 
kernel 0 ( 91 , 92 )^ 12 )^ over the spatial coordinates ri and r 2 as 

^ JJ n(ri)0(9i,92,ri2)n(r2)dridr2, (2) 

where ri 2 = |ri — r 2 l, 91 = 9 o(ri), 92 = 9 o(r 2 ), and 9 o(r) is a function of the charge density n(r) and 
its gradient |Vn(r)|. Roman-Perez and Soler (RPS)^ proposed an efficient method for evaluating 
E“*, as the direct computation of the double integral in Eq. ([2]) is extremely time-consuming. In 
the RPS scheme, the vdW kernel 0 is expanded in terms of the interpolating function PaiQ), which 






4 


satisfies PaiQp) = so that cp is a fixed value at a given point, allowing the use of the fast 
Fourier transform (FFT) in the evaluation of Fi”*. However, (p has a divergence when Qa-iQp —^ 0, 
which makes the interpolation near the origin difficult, and thus RPS replaced (p with the soft 
form, and an LDA-like approximation was employed for near the origin. Alternatively, Wu 
and Gygi (WG)^ proposed a simplified implementation to avoid the divergence in (p as follows. 
The nonlocal interaction kernel, multiplied by qi and q 2 , is expanded as 


QiQ2(p{qi,q2, ru) ^ qaPa{qi)qi3Pp{q2)(papiri2), (3) 

so that the divergence at qa^qp —^ 0 can be avoided. By introducing the function = 

gQ,n(r)p„((/o(r))/qo(r)j the nonlocal correlation is calculated in the reciprocal space as 

= jj Va{Yi)'nj3{v2)(pai3{ri2)dridr2 

ajS 

= (4) 

G a/3 

where G is the reciprocal lattice vector G = |G|, 3/a(G) and (pap{G) are the Fourier components 
of ? 7 a(r) and (paji{'''i 2 ), respectively, and 11 is the volume of a unit cell. The Fourier component 
(pa/siG) is calculated as 

<PMa) = psFaJ^), (5) 

\HOi / 

where 


i/ 3 (fc) ~ J siii{ku)du. 


( 6 ) 


Note that the G = 0 term in Eq. (|4]) vanishes when a = (5. We discuss the behaviors of the (p 
function in Appendix [Bj 

The nonlocal correlation potential within the White-Bird algorithm^^ is given by 

6EP:^ 


<(r*) = 


where 


5n{Yi) 

N dEf 

H dn[Ti) 

E 


11 


r,; 


,dr]a{ri) , ^ diqai^j) dVn{rj) 


dn{Ti 


+ ^Ua{rj) 


9Vn(rj) dn{ri) 




Pd 


/3G 


with N being the number of FFT grids. 


( 7 ) 


( 8 ) 
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B. Pressure tensor 


To optimize the lattice parameter, it is necessary to calculate the pressure tensor, and there is 
an additional contribution from the nonlocal correlation in vdW-DF.— >2^ In our formulation, an 
additional contribution to the pressure tensor (/c, £ = x,y, z) is calculated as 


(He 


n 

1 


dh 


km 


)rn£ 


= - 




dn{rj) dhkm 


<m£ 


+ 




dVa{rj) {Vn{Tj))k{Vn{rj))i 


N 

1 

iv 


EEE'^“(''. 

m j a 

Vn(rj) 


5|Vn(rj)| |Vn(rj)| 
dVairj) 


|Vn(r,-)| 


E 


^ d\Vn{rj)\ 

{h^UeiiG)E 




dh 


km 


dcj)^f,{G) GkGe 


a/? G III 


(9) 


( 10 ) 


where the matrix h = (ai,a 2 ,a 3 ) is a set of primitive lattice vectors a^ (fc = 1,2,3) and P is the 
inverse of h. 


C. Spin-polarization-dependent gradient correction of the local correlation 

In the original vdW-DF, LDA was employed for the short-range local correlation energy, to avoid 
possible double counting of the contribution from |Vn| contained in E^^.— In a spin-polarized sys¬ 
tem, spin-polarization-dependent gradient correction should be included in the correlation energy 
and potential, but such a contribution is missing in the nonlocal correlation, as it is mostly for¬ 
mulated for the spin-unpolarized system. In our previous work22^ we propose a simple scheme to 
include such spin-polarization-dependent gradient contribution as follows. We define the local part 
of the correlation functional as 


Fi°= = F;f°A[nt,nJ + AF;c[n,C], 


( 11 ) 
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where 


AE,[nX] = E^^^[n^,rn] - E^^^[n/ 2 ,n/ 2 ] 

= j n{H{rsX,t) - H{rs,0,t)} dr, (12) 

H, Vs, C) and t are the gradient contribution, Seitz radius (n = S/dvrr^), relative spin polarization, 
and dimensionless density gradient proportional to |Vn|, respectively. We use the functional form 
for H proposed by Perdew, Burke, and Ernzerhof (PBE).— The present functional is reduced to 
the original one in the absence of spin polarization. 

D. Technicalities 

We have implemented vdW-DE in our in-house DET code,— which uses a plane-wave basis set 
and ultrasoft pseudopotentials.— >2^ In the construction of the pseudopotentials, we neglected the 
nonlocal correlation (E^^) and employed the semilocal exchange and correlation functionals, which 
are consistent with the solid-state calculation. For instance, the pseudopotential constructed using 
the reht Perdew Wang (PW86R)^>^ exchange and the LDA correlation was used in vdW-DF2 
calculation. The use of the pseudopotential generated with the semilocal exchange correlation 
functional has been justified.—*^ We also performed the calculations using the pseudopotential 
generated with the PBE functional and found that there is an insignificant difference in binding 
energy (~ 1 meV/atom). In the present vdW-DF calculations, we used relatively large energy 
cutoffs for wave functions and electron density in the plane wave expansions, compared with those 
used in LDA or GGA, to achieve the convergence of the pressure tensor, because the vdW forces are 
very weak, and the potential energy surface is very sluggish. Unless otherwise noted, we optimized 
the lattice parameters and associated internal structural parameters to converge within a pressure 
of 0.05 GPa. The Brillouin zone integration was performed using the set of Monkhorst-Pack (MP) 
special fc-points.— In this work, we employed the Perdew-Zunger exchange correlation for LDA— 
and the PBE functional for GGA—. 

For the vdW-DF calculations, we employed the following functionals: the original vdW-DF with 
the revPBE^ exchange and LDA correlation, the second version of the vdW-DF (vdW-DF2), which 
uses PW86R— exchange, vdW-DF paired with Cooper’s (C09}^ exchange (vdW-DF*^^®’’^), and 
vdW-DF2 paired with C09 exchange (vdW-DF2*^*^®’’^). We also implemented the revised Vydrov- 
Van Voorhis functional (rVVlO) ^^'^^ with the WG^ implementation (see Appendix |A|), which uses 
the PW86R exchange, PBE correlation, and the nonlocal correlation. The functional is designed 
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TABLE I. Optimized lattice parameter (ofcc), equilibrium volume (Vq), and cohesive energy (AEcoh) for 
solid argon in the fee structure, obtained using different density functionals. The CCSD(T), ACFDT-RPA 


using PBE orbitals, and experimental values are also listed for comparison. 


Functional 

f^fee 

(A) 

Vo AFcoh 

(A^/atom) (meV/atom) 

LDA 

4.97 

30.61 

135 

PBE 

6.04 

55.03 

14 

vdW-DF 

5.53 

42.34 

142 

vdW-DF2 

5.29 

37.05 

116 

vdW-DF^o^’' 

5.34 

38.06 

105 

vdW-DF2C09’^ 

6.32 

63.06 

21 

rVVlO 

5.20 

35.20 

99 

rVV10“ 

5.17 

34.55 

117 

CCSD(T)^ 

5.251 

36.196 

87.9 

ACFDT-RPA'= 

5.3 

37.22 

83 

Expt. 

5.311'^ 

37.451^ 

80.R 


Ref. 


45. Ref. 


46. Ref. 


47. ^ Ref. 


48. Ref. 


49 


to give accurate Cg coefficients and to vanish in the uniform electron gas limit. In the case of 
antiferromagnetic a-02, we employed the spin-polarization-dependent gradient correction (GC) as 
described previously, and “-GC” is appended to the functional name when GC is considered^i. 
Note that GC is not necessary for rVVlO, because the PBE-GGA correlation is used for the local 
correlation. 


To evaluate (f> in the vicinity of = 0, we use a linear mesh up to g = qm with grid points, 
and the logarithmic mesh is used from Qm to Qc with grid points, where Qc is the cutoff for the 
g'-mesh. We use Qc = 8 a.u., qm = 0.01 a.u.,A^'“ = 3, and = 28, except for rVVlO, in which 
qc = 8 a.u., and qm = 0.005 a.u. are used. We confirmed that the cutoff qc is sufficiently large, 
by monitoring the distribution of go the nonlocal correlation energy as a function of go (see 
Appendix ICl) . 
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TABLE 11. Optimized lattice parameters (ohcp and Chcp), equilibrium volume (Vq), and cohesive energy 


(AEcon) for hep argon. 


Functional 

Q-hep Chep/^hep 

(A) 

To 

(A^/atom) 

AEcoh 

(meV/atom) 

LDA 

3.51 

1.633 

30.45 

135 

PBE 

4.16 

1.633 

50.92 

13 

vdW-DF 

3.89 

1.633 

41.63 

142 

vdW-DF2 

3.74 

1.636 

37.08 

116 

vdW-DFC°9’^ 

3.80 

1.632 

38.75 

105 

vdW-DF2C09^ 

4.19 

1.632 

52.81 

20 

rVVlO 

3.66 

1.633 

34.70 

99 

Expt.“ 

3.8 

1.63 

39 



a 


Ref. 
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III. RESULTS 

A. Fee and hep argons 

We begin with solid argon in the face-centered-cubic (fee) and hexagonal-elose-paeked (hep) 
struetures, as a typieal vdW bonded solid. It is known that at ambient pressure, argon in the fee 
strueture is the most stable, whereas that in the hep strueture is metastable. For both phases, we 
optimized lattiee parameters using several vdW-DFs as well as LDA and GGA. The energy eutoffs 
of 60 and 500 Ry were used for wave funetion and eleetron density, respeetively. A 14 x 14 x 14 
(8x8x8) MP speeial fc-point set was used for the fee (hep) phase. The optimized lattiee eonstant 
and eohesive energy for fee argon are summarized in Table HI along with GGSD(T) and experimental 
values. The lattiee eonstant obtained using LDA (PBE) is underestimated (overestimated),whereas 
those with vdW-DFs are improved exeept for vdW-DF2^*^®^, and most of the vdW-DFs improve 
eohesive energy. Our results for fee argon are in good agreement with the reeent results by Tran 
and flutter.— 

The results for hep argon are summarized in Table |TI1 The differenee between the eohesive 
energies for the fee and hep phases is very small (at most 2 meV/atom) and these phases are almost 
degenerated. The lattiee parameters for the hep strueture satisfy the ideal ratio Ofee ~ \/2ahcp) 
Chep/flhep = 1-633, where Ofee, Ohep, and Chep are the lattiee eonstant for the fee strueture, the 
in-plane lattiee eonstant for the hep strueture, and the out-of-plane-lattiee eonstant for the hep 
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structure, respectively. The deviations of a^cc from the ideal value are 0.15 A for PBE, 0.36 
A for vdW-DF2^^®^, and 0.03 A for other functionals. Our result indicates that these two phases 
are almost identical and cannot reproduce the experimentally observed stable fee phase at ambient 
pressure. This apparent discrepancy may be solved by introducing an entropic effect from phonons. 
Indeed, Ishikawa et al— performed lattice dynamics calculations and successfully predicted a stable 
fee argon by taking into account the entropic contributions. 


B. Graphite 

Graphite consists of the two-dimensional carbon allotrope, graphene, and graphene layers bound 
in the out-of-plane direction with a vdW interaction. Graphite has been studied extensively both 
experimentally and theoretically, and is often used to assess the accuracy of a new method for weak 
interactions, as benchmark calculations with a highly accurate electronic structure method such as 
QMG^ and ACFDT-RPA.^ 

In our calculation, the Brillouin zone was sampled using an 8 x 8 x 4 MP /c-point set, and plane 
wave cutoffs of 80 and 480 Ry were used for wave functions and electron density, respectively. We 
fully optimized the lattice parameters according to the calculated internal pressure, and binding 
energy was determined from the difference in total energy at the equilibrium and that at an 
interlayer distance larger than 13 A. In Table Hill optimized lattice parameters and binding energies 
obtained using different exchange-correlation functionals are summarized. 

In accordance with the literature, LDA gives lattice parameters, which are in good agreement 
with the experimental results, whereas GGA yields a much larger out-of-plane lattice parameter 
(c) and negligible binding energy, suggesting that the latter cannot predict the binding of graphite. 
The binding energy obtained with LDA is smaller than that obtained from the experiment and by 
a highly accurate theoretical method. In general, our vdW-DF results are in good agreement with 
those of previous studies, that employed different vdW-DFs.— vdW-DF and vdW-DF2 
predict the binding energy (AFb), in good agreement with the experimental results. They also im¬ 
prove the description of structural properties, but the equilibrium volumes (Vq) are overestimated. 
Vo is improved using vdW-DF^^®^ and rVVlO, but they overestimate AFb. We note that our AFb 
obtained using rVVlO is almost twice as large as that reported by the original authors.— At present, 
the origin of the discrepancy is yet to be clarified, but similar larger values (~70 meV/atom) were 
obtained using different implementationsSS.^^ of rVVlO and different potentials.— On the other 
hand, vdW-DF2^^®^ predicts that lattice parameters and AFb are in good agreement with the 
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TABLE III. In-plane lattice constant (a), out-of-plane lattice constant (c), equilibrium volume (Vb), and the 
binding energy AEb for graphite. The numbers in parentheses were obtained using the experimental value. 


Functional 

a 

(A) 

c Vo AFb 

(A) (A^/atom) (meV/atom) 

LDA 

2.45 

6.69 

8.71 

24 

PBE 

2.47 

8.76 

11.61 

2 

vdW-DF 

2.49 

7.18 

9.61 

55 

vdW-DF2 

2.48 

7.05 

9.41 

53 

vdW-DFC09’^ 

2.47 

6.50 

8.60 

76 

vdW-DF2C09’^ 

2.47 

6.58 

8.71 

56 

rVVlO 

2.48 

6.76 

8.93 

68 

rVV10“ 

2.46 

6.72 

8.80 

39 

VVIO^ 

(2.46) 

6.777 

(8.88) 

71 

vdW-DF-cx= 

2.46 

6.43 

8.54 

66 

QMC^' 

(2.46) 

6.852 

(8.98) 

56±5 

ACFDT-RPA® 

(2.46) 

6.68 

(8.75) 

48 

Expt. 

2.46^ 

6.70^ 

8.80-^ 

52 ±59 


“ Ref. 


26. Ref. 


. Q 


Ref. 


27. Ref. 


52 


Ref. 


53. f Ref. 


55. 9 Ref. 


56 



FIG. 1. (Color online) Trigonal structure of selenium. 


experimental results, in line with a previous study.— 
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TABLE IV. Optimized lattice parameters (a and c), equilibrium volume (Vq), internal parameter (x), shortest 
Se-Se distance in a chain (^i), shortest Se-Se distance in different chains (£ 2 ), and binding energy (AEh) for 
trigonal selenium. 


Functional 

a 

(A) 

c 

(A) 

Vo 

(A^/atom) 

X 

£1 

(A) 

£2 

(A) 

AFb 

(meV / atom) 

LDA 

3.89 

5.04 

21.94 

0.2540 

2.40 

3.06 

395 

PBE 

4.51 

4.97 

29.16 

0.2148 

2.36 

3.58 

55 

vdW-DF 

4.70 

5.04 

32.15 

0.2091 

2.39 

3.74 

189 

vdW-DF2 

4.57 

5.12 

30.81 

0.2177 

2.42 

3.62 

228 

vdW-DFC°9’^ 

3.96 

5.11 

23.08 

0.2521 

2.43 

3.11 

429 

vdW-DF2C09^ 

3.99 

5.10 

23.42 

0.2497 

2.42 

3.14 

330 

rVVlO 

4.24 

5.11 

26.54 

0.2348 

2.42 

3.35 

320 

Expt. 

4.366“ 

4.954“ 

27.261“ 

0.225'' 

2.373“ 

3.436“ 



a 


Ref. 


64. 


Ref. 


65 


C. Trigonal selenium 

Trigonal selenium is formed by a bundle of one-dimensional covalently bonded atomic chiral 
chains. As already shown in previous works,— the gradient correction to LDA improves the de¬ 
scription of selenium, but is still less satisfactory, presumably because of the lack of vdW interaction 
between the chiral chains in GGA. Bucko et ai— demonstrated using the semi-empirical dispersion 
correction that the structural parameters are significantly improved, and the agreement with the 
experimental results becomes much better. Here, we use vdW-DF to address the importance of 
vdW interaction. 

The trigonal selenium belongs either to the P3i21 or P3221 space group and has three atoms per 
unit cell (see Fig. [1]). The atomic position in a cell can be described using the internal parameter 
X, which scales the distance from the screw axis to the atomic position. The structure of the 
trigonal selenium is also characterized by the shortest Se-Se distance in the same chain (£ 1 ) and 
the shortest Se-Se distance in different chains (£ 2 )- We used an 8 x 8 x 4 MP fc-point set and 
plane wave cutoffs of 60 and 500 Ry for wave functions and electron density, respectively. Binding 
energy was calculated from the difference between the total energies of a solid in equilibrium and 
the isolated chain. 

Calculated structural parameters and binding energy are summarized in Table IIVI In general. 
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FIG. 2. (Color online) Cubic Pai structure of dry ice (solid C02). 


the lattice constant along the chiral chain (c) is overestimated by all the functionals used in this 
study, whereas the accuracy for the lattice constant a varies. LDA significantly underestimates a, 
whereas PBE significantly improves the description of a, in good agreement with previous studies. 
Both vdW-DF and vdW-DF2 overestimate a and c, whereas the use of the COD exchange (vdW- 
P)pC09x vdW-DF2^^®^) leads to the underestimation of a. On the other hand, although c is 
slightly overestimated, rVVlO provides a balanced description of both the lattice constants. It is 
found from our results that there is a correlation between calculated a and c: if a is overestimated, 
the error in c tends to be smaller. Thus, to obtain accurate structural parameters for the trigonal 
selenium, it is very important to describe both covalent and vdW interactions accurately, suggesting 
that this material can be a good benchmark system for a vdW inclusive functional. 

D. Dry ice 

The solid form of carbon dioxide is called dry ice. Carbon dioxide has no dipole moment, and 
thus attractive vdW forces play an important role in the condensation. We chose this material as 
a representative application of our vdW-DF implementation to non magnetic molecular crystals. 

Crystalline dry ice has a cubic symmetry with the space group of Pa3 (see Fig. [2] for the 
structure). Plane wave cutoffs of 160 and 960 Ry were used for wave functions and electron 
densities, respectively. An 8x8x8 MP /c-point set was used for Brillouin zone sampling. 

Optimized structural parameters and binding energy for the dry ice are given in Table IVj 
along with the MP2^ and experimental (I50K)^ results. Our PBF equilibrium volume is in good 
agreement with the theoretical value obtained by Bonev et ai— (52.9 A^/C02) using the same 
functional. On the other hand, both the equilibrium volume and binding energy obtained using 
vdW-DFs are in better agreement with the highly accurate MP2 and experimental (I50K) results, 
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TABLE V. Equilibrium lattice parameter (a), volume (Vq), C-0 bond length (£b), and binding energy (AEb) 


for dry ice. 


Functional 

a 

(A) 

Po 

(A3/CO2) 

4 

(A) 

AEb 

(meV/COz) 

LDA 

5.28 

36.88 

1.165 

370 

PBE 

6.04 

55.11 

1.175 

103 

vdW-DF 

5.77 

47.94 

1.180 

364 

vdW-DF2 

5.61 

44.04 

1.178 

346 

vdW-DF^°®’" 

5.53 

42.17 

1.176 

339 

vdW-DF2C09’' 

5.79 

48.42 

1.176 

155 

rVVlO 

5.51 

41.72 

1.179 

328 

MP2“ 

5.46 

40.69 

1.17 

290 

Expt. 

5.62 '' 

44.38'’ 

1.155'’ 

288'’ 


a 


Ref. 


69. 



suggesting the improvement over LDA and PBE. The maximum deviation of the equilibrium volume 
obtained using vdW-DF is 9.1% compared with that obtained in the low-temperature experiment. 
The above results suggest that PBE overestimates the equilibrium volume because of the lack of 
dispersion forces, and a more accurate description of the structure and energetics of dry ice is made 
possible by considering the dispersion forces with vdW-DE. Regarding the severe underestimation 
of the binding energy using vdW-DE2^^®^, it has been found^ from a systematic assessment using 
the S22 dataset that although it predicts reasonable intermolecular separation (distance), the 
functional severely underestimates the binding energy, because the 009 exchange is too repulsive 
at a relatively large density gradient relevant to the intermolecular region. This result implies that 
vdW-DF2‘^^®^ inaccurately describes the intermolecular vdW interaction. 


E. Solid oxygen in the a phase 

a -02 has an antiferromagnetic ground state with the crystal structure of the C2lm space group, 
as shown in Fig. [3l There are four lattice parameters, a, b, c, and /3, and the internal parameters of 
ih (bond length of O 2 molecule) and 9 (tilted angle of molecular axis). The molecular axis is almost 
perpendicular to the a6-plane and slightly tilted within the ac-plane. The angle 9 was reported to 
be ~3° in the experiment.—. In the calculation, we used plane wave cutoffs of 160 and 960 Ry for 
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FIG. 3. (Color online) Crystal structure of solid oxygen in a phase. The arrows on molecules represent the 
magnetic moment. 


TABLE VI. Optimized lattice parameters (a, b, c, and /3), nearest-neighbor distance {imoi = \/f?T^/2), 
equilibrium volume (Vo), bond length (ih)-, magnetic moment {Ma) on oxygen atom, and binding energy of 


molecule (AEb). Experimental values are shown for comparison. 


Functional 

a 

(A) 

b 

(A) 

c 

(A) 

/3 

(deg) 

^mol 

(A) 

Do 

(A3) 

h 

(A) 

Ma 

(he) 

AEb 

(meV) 

LDA 

3.29 

3.28 

4.05 

113.9 

2.32 

19.93 

1.202 

0.30 

552 

PBE 

4.59 

3.93 

5.05 

122.1 

3.02 

38.54 

1.218 

0.66 

41 

vdW-DF 

4.68 

3.68 

4.70 

125.2 

2.98 

33.05 

1.231 

0.66 

221 

vdW-DF-GC 

4.94 

3.57 

4.91 

128.4 

3.05 

33.85 

1.231 

0.66 

213 

vdW-DF2 

3.83 

3.85 

4.29 

118.6 

2.72 

27.75 

1.235 

0.60 

225 

vdW-DF2-GC 

3.91 

3.88 

4.30 

119.0 

2.76 

28.56 

1.235 

0.62 

209 

vdW-DF^o^’' 

3.52 

3.46 

4.15 

114.8 

2.47 

22.90 

1.215 

0.47 

285 

vdW-DFC09^-GC 

3.59 

3.47 

4.17 

115.4 

2.50 

23.51 

1.215 

0.50 

255 

vdW-DF2C09x 

3.66 

3.56 

4.34 

115.1 

2.55 

25.51 

1.216 

0.52 

95 

vdW-DF2C09x_Gc 

3.79 

3.62 

4.36 

115.8 

2.62 

26.92 

1.217 

0.57 

71 

rVVlO 

3.71 

3.62 

4.18 

117.2 

2.59 

24.94 

1.225 

0.56 

240 

Expt.“ 

5.403 

3.429 

5.086 

132.3 

3.200 

34.85 

1.29 




“ Ref. 


4. 


wave functions and electron density, respectively, and an 8 x 8 x 8 MP /c-point set was used for the 
Brillouin zone sampling. During the structural optimization, the molecular axis was fixed in the 
direction perpendicular to the a6-plane, but the residual forces on the atom were small, typically 
less than 0.001 eV/A. 

In Table IVIl we report the structural parameters and binding energies obtained using different 
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TABLE VII. Magnetic energy for a-02 


Functional 

(meV) 

PBE 

95 

vdW-DF 

118 

vdW-DF-GC 

87 

vdW-DF2 

304 

vdW-DF2-GC 

251 


functionals. As expected, LDA severely underestimates the eqnilibrinm volume by 42.8%, whereas 
PBE overestimates it, but the error is marginal (10.6%). However, this unexpectedly small error 
is dne to the cancellation of the errors in a (—15.0%) and b (14.6%). The error in c is surpris¬ 
ingly small, but the trend in the lattice parameters estimated using PBE is not systematic. On 
the other hand, all the vdW-DFs underestimate a and c, whereas b is overestimated, and as a 
result, the equilibrium volumes are consistently underestimated. The angle /3 is also consistently 
underestimated. Among vdW-DFs, the original vdW-DF proposed by Dion et al. predicts the 
most accurate structural parameters, and inclnsion of GC further improves the structural parame¬ 
ters, which are in good agreement with the experimental results, suggesting the importance of the 
spin-polarization-dependent gradient correction of the local correlation. 

The nearest-neighbor distance between O 2 molecules (.^moi = + ^^/2) is a very impor¬ 

tant quantity in 0 - 02 , as it strongly correlates with the antiferromagnetic interaction between O 2 
molecules. We found that is nnderestimated by LDA and PBE, but is increased using vdW- 
DF-GC to a distance of 3.05 A which is close to the experimental value (3.20 A). We also note that 
there is a correlation between /3 and c estimated using vdW-DF, i.e., the smaller the c, the smaller 
the /3. As discussed in previons works,— >2^ the distance between magnetic molecules is determined 
by a subtle balance between magnetic and vdW interactions. In the solid state, not only the bal¬ 
ance between magnetic and vdW interactions within the a6-plane, but also that between a6-planes 
(out-of-plane direction) plays a decisive role: a complicated interplay among the antiferromagnetic, 
ferromagnetic, and vdW interactions in three dimensions determines the structure of a-02. The 
molecules in the a6-plane can be considered to form a triangular configuration. In the experimental 
structure, this configuration is very similar to an equilateral triangle, where the antiferromagnetic 
interaction between the neighboring molecnles could destabilize the magnetic structure or distort 
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the triangle to reduce the magnetic frustration. On the basis of this consideration and the knowl¬ 
edge acquired in a previous stud}*^^, the underestimation of the nearest-neighbor distance (.^moi) 
using vdW-DF may be ascribed to the overestimation of the antiferromagnetic interaction between 
molecules with respect to the ferromagnetic one. The structural properties depend strongly on the 
exchange energy functional used, as shown in Table IVIl Such a behavior has also been shown in 
the previous study on the pair of oxygen molecules^^. The revPBE exchange functional in the orig¬ 
inal vdW-DF gives a more repulsive intermolecular potential than the other exchange functionals. 
The other vdW-DFs consisting of different exchange functionals (C09 and PW86R) underestimate 
^moi) indicating that these exchange functionals enhance the overstabilization of the antiferromag¬ 
netic state. Thus, to improve the description of a-02, it is necessary to develop more accurate 
exchange and correlation functionals that predict a balanced description of antiferromagnetic and 
ferromagnetic interactions. 

The magnetic interaction J is proportional to —t'^/AE, where t and AE are the transfer integral 
and the energy gap between the spin-up and spin-down states immediately above and below the 
Fermi level, respectively. Because GGA is known to underestimate the energy gap for insulating 
and molecular systems, the antiferromagnetic states are over stabilized as a result of the underesti¬ 
mation of AE (overestimation of |t| as well). To capture a trend in the functionals on the magnetic 
interaction, we calculated magnetic energy, defined as the energy difference between the ferromag¬ 
netic (F) and antiferromagnetic(AF) states at the same AF crystal structure = E^ — E^^). 

The results obtained using PBE, vdW-DE, and vdW-DE2 (with and without GC) are summarized 
in Table IVIIl It is found that vdW-DF, which predicts more accurate lattice parameters (larger 
^moi) than vdW-DF2, yields a smaller suggesting that the intermolecular distance plays 

an important role in the magnetic interaction. Inclusion of GC leads to the enlargement of the 
nearest O 2 distance (.^moi) and reduction of AE™^®, because GC destabilizes the antiferromagnetic 
states, which works as a repulsive force between the nearest O 2 molecules22. 


IV. SUMMARY 

In this work, we have studied solid oxygen in the a-phase with the antiferromagnetic configura¬ 
tion, using several variants of vdW-DE and the correction scheme for magnetic systems proposed 
in our previous study. To test the vdW-DF energy, force, and stress implemented in our program 
code, we performed the calculations for solid argon in the fee and hep phases, graphite, trigonal 
selenium, and dry ice, in which the vdW interaction is considered to be important. We have found 
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that vdW-DF shows an overall improvement in the structure and energy of these materials over 
LDA and GGA, in good agreement with previous studies. In the case of antiferromagnetic solid 
oxygen, we have found that vdW-DFs consistently underestimate the lattice parameters, and the 
original vdW-DF proposed by Dion et al. with our spin-polarization-dependent gradient correction 
scheme (vdW-DF-GC) provides the most accurate structural parameters among other functionals. 
We point out the competing nature of ferromagnetic, antiferromagnetic, and vdW interactions in a 
solid oxygen, and for a more accurate prediction of the structural properties of solid oxygen, a bal¬ 
anced description of these interactions is decisively important. Thus, we suggest that solid oxygen 
can be a critical test material for an electronic structure method, which can account for magnetic 
and vdW interactions accurately. Applications of the present vdW-DF scheme for the magnetic 
systems to the molecular magnet or metalorganic molecular systems are highly anticipated, which 
promote research on molecular spintronics. 
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Appendix A: rVVlO with the Wu-Gygi implementation 


The rVVlO functional has been proposed^^ for implementing the VVIO functional^ with the 
efficient RPS algorithm.— We note that, very recently, the original VVIO has been implemented 
using the WG scheme by Corsetti et ali^ Here, we briefly describe our implementation of rVVlO 
with the algorithm proposed by WG.— The rVVlO nonlocal correlation functional is given by 




^(g) ^rVVlO/- - ^ ^(g) j j 

-^ (gEg2,n2) ,.,/o, ^ Qridr2, 


(Al) 


2JJ A:3/2(ri)^ ""'A:3/2(r2) 

where ^ 2 ) ^ 2 ) is the rVVlO nonlocal correlation kernel,— qi = q{ri), q 2 = 9 (i' 2 )) and 

k{n{^ = 37r6(n/97r)^/®/2. qi is defined as q{ri) = a;o(n(rj), |Vn(ri)l)//c(n(rj)) with uq defined in 
Ref. [23|. The parameter b is an adjustable parameter, which is determined by fitting on a training 
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FIG. 4. (Color online) Fourier component of the kernel function (j)ap{G) for several sets of qa and qp. 
Diagonal (upper panel) and off-diagonal (lower panel) components are shown. 


set (see Refs. 


23 and 


261 for details). By expanding the nonlocal correlation kernel in the same 


manner as Eq. the nonlocal correlation energy is written as 


=§EE^«(G)77;3(G)</-™°(G)> (A2) 

G q:/3 

where ? 7 a(G) is the Fourier component of f)a(r), which is defined by fia{r) = qan{r)pa{q{r))/{r)q{r), 
and the Fourier transform of the nonlocal correlation kernel jg gjygj^ j^y 

(A3) 

FqY^^°(/c) f ( u, du. (A4) 

^ JO \ \ Qa / 

The nonlocal correlation potential is calculated similarly to Eq. (j7|). 


Appendix B: Fourier component of the kernel function cj)ap{G) 

In the present implementation, the Fourier components of the kernel function (j)ap{G) are cal¬ 
culated from the tabulated data for 4>{di, ^ 2 ) {di(t>{di,d 2 ) in the case of WG implementation) using 
Eq. dH). The typical functions of 4’a/3iG) are presented in Fig. 01 For the small q^s, (paa (G) 
deviates slightly from zero at G = 0 in the present computation, which may cause a numerical 
error. However, the error in the absolute value of F))* associated with the deviation is typically 
~0.5 meV in fee argon and negligible in the estimation of cohesive energy. 
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FIG. 5. (Color online) Distributions of go function (solid curves) and ec*(a) (dotted curves) for fee argon 
calculated using vdW-DF (upper panel) and vdW-DF2 (lower panel). The red open and blue filled symbols 
are for the lattice parameters (a) of 5.29 and 10.58 A, respectively. 


Appendix C: Distribution of go function 

To expand Eq. ([3]), an appropriate cutoff qc and mesh for the go function (qa) are needed. To 
determine qc and g^, we investigated the distribution of the go function by sampling the values 
on each EFT grid point. The distribution Tl(go) was calculated as a histogram by counting the 
number of samples while satisfying g^ < go(r) < go+i- In Fig. [5l Tl(go)’s for the vdW-DF and 
vdW-DF2 are plotted for different lattice parameters. Furthermore, we analyzed the contribution 
of the nonlocal correlation energy from each g^ value, by defining 

a 

^c'(«) = § E ( \Vc.{G)\^cf>caiG) + 2 ^ Re«(G)r?^(G))0„^(G) ) , (Cl) 

G \ /3«a) / 

where ^“^(a) is a positive value. e“*(a)’s for fee argon calculated using vdW-DF and vdW-DF2 
are shown in Fig. [SJ There is almost no contribution of the nonlocal correlation energy at the 
smallest g^, although the samples of D{qa) exist there. Note that the samples at the smallest g^ 
come from diluted electron densities. In the case of a = 5.29 A, e“*(a)’s have a peak at g = 2.4 
and are damped rapidly as g^ increases. In the case of a = 10.58 A, the shrunken electron cloud 
around atomic cores causes the distributions of D{qa) at larger gQ,’s(> 5). However, there is no 
energy contribution from the term. This is presumably because the vdW interaction is both 
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canceled out in the atom and fully damped at the atomic distances. From the above result, it is 
confirmed that the qa mesh defined in the present work covers the energy range in the system. 
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